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In this study, we present an integrative systematic revision of the spider gecko, Agamura senso lato, in lran. We sampled 56 
geckos of this complex from its distributional range in Iran and western Pakistan and sequenced these for two mitochondrial 
markers, cytochrome b and 12S ribosomal RNA, and one nuclear marker, melano-cortin 1 receptor. We combined our 
molecular data with species distribution modelling and morphological examinations to clarify Agamura persica systematics 
and biogeography. Due to a lack of published data, we used only our data to investigate the spatial and temporal origin of 
spider geckos within a complete geographic and phylogenetic context. The phylogenetic analyses confirm the monophyly of 
Agamura. Among spider geckos, Rhinogekko diverged around the early-mid Miocene (17 Mya) from the Lut Block, and then 
Cyrtopodion diverged from the Agamura clade about 15 Mya in the mid-Miocene as a result of the uplifting of the Zagros 
Mountains. Subsequent radiation across the Iranian Plateau took place during the mid-Pliocene. Agamura kermanensis 
exhibits deep divergence from two other species of Agamura (A. persica and A. cruralis), whereas no geographical substructure 
was observed on the Iranian Plateau for A. persica and A. cruralis. Our findings reveal that diversification is consistent with 
a biogeographical model explained by different dispersal waves and vicariant events on the Iranian Plateau during the last 
18 Mya. The divergence times between clades are compatible with orogenic events in southern Iran that resulted from the 
collision with Arabia. According to the genetic differentiation of both mtDNA genes (12S and cytochrome b), the systematic 
status of A. cruralis is confirmed, the new clade was distinguished from the genus Agamura, monophyly of Rhinogekko was 
confirmed and the allocation of Cyrtopodion gastrophole to the Cyrtopodion clade was confirmed. 
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INTRODUCTION in the region such as mesic, alpine and xeric, which are 
restricted by mountain chains (Farahmand & Nazari, 
he genus Agamura has awide distribution rangeinthe 2015). The formation of the plateau started ca.40 million 
Iranian Plateau and occupies different habitat types Years ago (Mya) as a result of the merger/collision among 
(Anderson, 1999). However, due to the different habitat | Arabian, Eurasian and Indian plates (Mouthereau, 2011). 
types, a high degree of geographic variation and genetic The formation and the geological history of the Iranian 
divergence exist. The Iranian Plateau is located in south- Plateau has had a fundamental influence in shaping 
western Asia and is comprised of various ecoregions __ the distribution patterns of the reptilian biota and the 
(e.g., both mountains and dry deserts), surrounded by _— SPeciation processes by both dispersal and vicariance 
several mountain ranges such as the Zagros in the west, events (Macey et al., 2000). 
Alborz in the north, Kopet Dagh in the north-east and Iran is an interesting region to study biogeographical 
Hindokush and Soleiman in the east and the south-east Patterns because it is located near the junction of 
(Macey et al., 1998). These mountain ranges create | the Eurasian, Arabian and Indian continental plates 
a rain shadow, preventing high precipitation on the (Macey et al., 1998). Collisions between these plates 
plateau, creating hot dry deserts in this central region created different mountain chains that affected the 
(Ahmadzadeh et al., 2012). The Iranian Plateau has a _—‘herpetofaunal biodiversity in the region (Van Hinsbergen 
high number of endemic lizards, especially within the arid et al., 2012; Anmadzadeh et al., 2012). In south-western 
clades of geckos (Smid et al., 2014). This highendemicity Asia, several phylogenetic studies have been conducted 
may have been facilitated by the various habitat types to gain insights on the biogeographical patterns and 
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et al., 2013; de Pous et al., 2015). One such example is 
the study of Smid et al. (2013) on the genus Hemidactylus 
and the following descriptions of several new species. 
Bauer et al. (2013) presented the phylogeny of naked- 
toed geckos suggesting that most speciation events in 
the western Palearctic took placed during the Miocene 
by vicariance due to different uplift events in the Iranian 
Plateau. There are more than 15 genera of the naked- 
toed gecko and one of them is Cyrtopodion that was later 
divided into several genera (Szczerbak & Golubev, 1996; 
Bauer et al., 2013). 

Within Gekkonidae, the spider gecko, Agamura, is a 
monotypic genus, phylogenetically close to Cyrtopodion, 
Rhinogekko, Bunopus and Crossobamon (Bauer et 
al., 2013), which ranges across Iran, Pakistan and 
Afghanistan. In Iran, it is known throughout most of the 
Iranian Plateau, east of the Zagros Mountains and south 
of the Alborz and Kopet Dagh mountain ranges, but is 
absent from the Kavir and Lut central deserts. Agamura 
occupies stony and rocky habitats, including hillsides and 
barren plains, with sparse shrubby vegetation (Anderson, 
1999, Sindaco & Jereméenko, 2008). To date, neither 
morphological investigations nor molecular studies on 
this monotypic genus have been carried out across its 
distribution range. Agamura has previously been used 
as an outgroup in other studies (Cervenka et al., 2010; 
Bauer et al., 2013; de Pous et al., 2015); however, these 
studies did not investigate the phylogenetic relationships 
within the genus, or its population structure and diversity. 
Agamura previously included four taxa (i.e., A. persica, A. 
gastropholis, A. misonnei and A. femoralis; (Szczerbak & 
Golubev, 1996), all of which were endemic to the Iranian 
Plateau. Later studies examined the genus and excluded 
the three latter species from Agamura: A. misonnei and 
A. femoralis were allocated to the genus Rhinogekko, 
and A. gastropholis to the genus Cyrtopodion (Anderson, 
1999; Krysko et al., 2007). Thus, currently the genus 
Agamura includes a single species, A. persica, which 
exhibits differences among populations in colour pattern 
and morphological characteristics (Szczerbak & Golubev, 
1996; Anderson, 1999). The subspecies A. p. persica and 
A. p. cruralis represent western and eastern populations, 
respectively (Anderson, 1999; Szczerbak & Golubev, 
1996). Differences between these two subspecies are 
apparent in the shape of the rostral scale, the number 
of the dorsal scales, and the sharper shape of the 
tubercles in the eastern populations (Anderson, 1999). 
Rhinogekko misonnei is another species that is endemic 
to Iran and only distributed around the Lut Desert and 
its phylogenetic status among other gekkonid species 
requires verification, as there has been no study on the 
species (Smid et al., 2014). 

In this study, we conduct an integrative revision of the 
monotypic genus Agamura senso lato in Iran. We provide 
the first comprehensive morphological and molecular 
analyses throughout its wide range in the Iranian 
Plateau with the aim of elucidating the evolutionary 
and biogeographic history of this enigmatic genus. For 
this purpose, we sampled different taxa throughout 
the Iranian Plateau including A. persica (different 
populations of the species from whole distribution range 


in Iranian Plateau) and its close relatives, Rhinogekko 
misonnei (providing the first genetic data for this genus), 
C. gastrophole, C. persepolense, and C. agamuroides. We 
used nuclear and mitochondrial markers to revise the 
systematics of Agamura. We used species distribution 
modelling and phylogeographic analyses to explore the 
phylogeographic structure of Agamura persica on the 
Iranian Plateau, and climate suitability that can affect the 
divergence between species in Iranian Plateau. 


MATERIALS AND METHODS 


We used two criteria to assess the species limits; first, 
the identification of lineages based on mitochondrial and 
nuclear markers and, second, the presence of diagnostic 
morphological characters. 


Sampling, DNA extraction and amplification 

A total of 56 individuals were collected from the Iranian 
Plateau during field trips from 2014 to 2016. We used 
two recently described species of the genus Cyrtopodion 
that is clearly situated in the main clade of Cyrtopodion 
(Nazarov et al., 2009; Nazarov et al., 2012). The dataset 
contains seven recognised species: Agamura persica, 
A. cruralis, A. kermanensis, Rhinogekko misonnei, 
Cyrtopodion gastrophole, C. persepolense, C. sistanense 
and C. agamuroides (Hosseinian Yousefkhani et al., 2018). 
Localities and coordinates for each sample are presented 
in the Supplementary Materials, Table $1. Specimens 
were deposited in the Sabzevar University Herpetological 
Collection (SUHC). Sequences of Hemidactylus turcicus 
were retrieved from GenBank and used as outgroup 
(Table S1). 

DNA was extracted from tissue samples using the salt 
method (Kabir et al., 2006). The quality of extracted DNA 
was measured using 1% agarose gels stained by 0.5 ul 
GreenViewer 6X and visualised under ultraviolet light. We 
amplified two mitochondrial genes 12S rRNA (12S) using 
the primers 12SL (5’-AAACTGGGAT TAGATACCCCACTAT- 
3’) and 12SH (5'-GAGGGTGACGGGCGGTGTGT-3’) (Kocher 
et al., 1989), and Cytochrome b (Cytb) with the primers 
L14724 (5'-GACCTGCGGTCCGAAAAACCA-3’) and H16064 
(5'-CTTTGGTTTACAAGAACAATGCTTTA-3’) (Burbrink et 
al., 2000) and L14919 (5’-AACCACCGTTGTTATTCAACT- 
3’) and Ei700r (5’-GGGGTGAAA GGGGATTTTRTC-3’) 
(Rastegar-Pouyani et al., 2010) and one nuclear gene 
Melano-cortin 1 receptor (MC1R) with the primers 
MC1R-F (5'-AGGCNGCCATYGTCAAGAACCGGAACC-3’) 
and MC1R-R (5'-CTCCGRAAGGCRTAAATGATGGGGTCCAC- 
3’) (Eskandani et al., 2010). 


Phylogenetic analyses and haplotype network 
construction 

We used Clustal W as implemented in Bioedit alignment 
editor v. 7.0 (Hall, 1999) to align sequences with default 
parameters. Protein coding sequences (Cytb and MC1R) 
were translated into amino acids with Mega v.6.0 (Tamura 
et al., 2013) and no stop codons were observed. Un- 
corrected genetic distances (p-distance) were calculated 
using Mega v.6.0 and ExcaliBAR (Aliabadian et al., 2014) 
for 12S and Cytb gene fragments independently. 


The best-fit models of nucleotide evolution were 
assessed using ModelTest 3.7 (Posada & Crandall, 1998) 
and the best fit models of evolution according to the 
Akaike Information Criterion (AIC) were: 12S-GTR+I+G; 
cyt b-TIM+1+G; mcir-TrN+I+G. 

The phylogenetic analyses were performed using 
Maximum Likelihood (ML) and Bayesian Inference (Bl) 
methods and for this purpose, all gene alignments were 
combined into a single alignment totalling 1670 bp (389 
bp of 12S; 625 bp of cytb; 656 bp of mcir). We also 
considered Hemidactylus as outgroup in the analyses 
(Pyron et al., 2013). We used 50 cytb sequences from 
Genbank to clarify the phylogenetic structure of the 
genera. 

Maximum Likelihood analyses were conducted 
with RaxML 7.4.2 (Stamatakis, 2006) as implemented 
in RaxmlGUI 1.3 (Silvestro & Michalak, 2012) with a 
GTR+I+G model. The analyses were run in heuristic 
search and the nodal support was obtained by bootstrap 
analysis with 1000 replicates (Felsenstein, 1985). 
MrBayes 3.2.1 (Ronquist et al., 2012) was used for the 
BI analyses and the best-fit models were specified above 
for concatenated dataset. The analyses were run for 
10’ generations with a sample frequency of every 1000 
generations. Some parameters, like the number of runs 
and the number of chains, were kept as default and a 
sufficient number of generations were evaluated by 
the log likelihood value (InL) and split frequency lower 
than 0.01. We conservatively discarded the first 25% of 
trees as burn-in (Condamine et al., 2015). To reconstruct 
the ancestral area Bayesian Binary MCMC (BBM; Ali et 
al., 2012), we employed Reconstruct Ancestral State in 
Phylogenies (RASP) (Yu et al., 2012) using all spider gecko 
sequences. MrBayes was used to prepare the input tree 
file. Six areas were designated based on zoogeographical 
regions for reconstruction as: Central Plateau, East Iran, 
South Iran, South-west Pakistan, South-west Iran, and 
Zagros Mountains. We chose these areas to identify the 
direction of dispersal within the Iranian Plateau spider 
geckos. 

Relationships among lineages and species were 
assessed with allele network of the mcir nuclear marker. 
The nuclear alignments were imported into TCS 1.21 
(Clement et al., 2000) using a parsimony method to 
obtain the haplotype network. 


Estimation of divergence time 

Because of the absence of internal calibration points 
for spider geckos and their relatives, we applied direct 
estimations obtained from other groups of lizards. The 
substitution rate of the same mitochondrial genes (12S 
and Cytb) that were calculated for three lizard families 
from the Canary Islands: Tarentola (Phyllodactylidae) 
(Carranza et al., 2000), Gallotia (Lacertidae) (Cox et 
al., 2010) and Chalcides (Scincidae) (Brown & Pestano, 
1998) were used to estimate the divergence time. These 
substitution rates have already been used for divergence 
time estimates for different taxa including Hemidactylus, 
Bunopus, Asaccus and etc. (Carranza & Arnold, 2012; 
Sindaco et al., 2012; Smid et al., 2013). BEAST 1.8 (Heled 
& Drummond, 2010) was used to estimate divergence 


Molecular phylogeny of spider geckos 


time among the spider geckos and the models and 
priors were applied as follows (otherwise by default): 
evolutionary models were set for each gene separately; 
random starting tree; clock models were set as lognormal 
relaxed clock with unlinked status; tree priors were set as 
coalescent and constant size. Finally, divergence times 
were assessed by the mean rate of molecular evolution 
for the ucld. priors for 12S (mean: 0.00755, stdev: 
0.00247) and Cytb (mean: 0.0228, stdev: 0.00806) gene 
fragments (Carranza & Arnold, 2012) independently. 


Species distribution modelling 

A total of 189 presence records from the examined species 
and clades were obtained from the literature, museum 
records and our direct filed surveys (Supplementary 
Materials, Table S2). Climatic layers were downloaded 
from the worldclim website (www.worldclim.org) in 30 
arc second (Hijmans et al., 2005) and extracted using 
ArcGIS 10.3 (ESRI) only for Iran (Table $3). 

Correlations between climatic layers were calculated 
usingENMTools1.3(Warrenetal.,2010)andthecorrelative 
layers (>0.7) were removed from the analyses. Maxent 
3.3.3e (Phillips et al., 2006) was employed to predict the 
potential distribution area using only presence records. 
The final set of variables with lower correlation than 0.7 
used for all species distribution modelling consisted of 
12 bioclimatic variables (Table 2). All models were run 
for 10 replicates under a crossvalidate model and 10000 
background points, with a convergence threshold of 
0.00001, and maximum number of iterations as 500. 
The model accuracy was evaluated using area under 
the curve (AUC) criterion that ranges between 0 and 1 
(Fielding & Bell, 1997). Model visualisations were done 
by ArcGIS 10.3 (ESRI) and we exported relevant maps as 
the species distribution prediction. 


Morphological analyses 

Populations of Agamura were examined using 27 
morphological characters (12 metric and 15 meristic 
characters; Supplementary Materials, Table S4) on the 
samples that were sequenced for phylogenetic analyses. 
Metric characters were measured using digital callipers 
(rounding to nearest 0.1 mm) and meristic characters 
were examined using an Olympus loupe. Operational 
taxonomic units (OTUs) were classified according to 
the cluster analysis and zoogeographic regions on the 
Iranian plateau (Anderson, 1999). Three OTUs were 
defined as western, eastern and southern clades. 
Analysis of variance (ANOVA) performed on the OTUs 
and the significant variables (P< 0.05) were excluded to 
run principal component analysis (PCA) and canonical 
variate analysis (CVA) and to visualise the morphological 
variation by these analyses. 


RESULTS 


Taxon sampling and sequence data 

Our dataset included mitochondrial fragments of 12S 
(389bp; V = 182; Pi = 136) and Cytb (625 bp; V = 297; 
Pi = 239) and a nuclear gene fragment MC1R (656 bp; 
V = 94; Pi = 63) totalling 1670 bp. Thirty-six unique 
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Figure 1. Sampling localities of spider geckos including all species. Numbers correspond to specimens listed in Table S1 and 
colours to specimens in Fig. 2 and 3. The black region in the right-top map indicates the sampling region in the Iranian Plateau. 


Table 1. Uncorrected genetic variation (p-distance) among different species of angular-toed geckos in the Iranian Plateau. 
Above the diagonal represents the variation in 12S and below the diagonal refers to Cytb diversity. 
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(1) A. cruralis 
(2) A. persica 


(3) A. kermanensis 
(4) C. agamuroides 
(5) C. persepolense 
(6) C. gastrophole 


(7) R. misonnei 
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13.2 
20.1 
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14.2 
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19.5 
19.1 


18.5 
20.1 
21.7 
21-1 


Table 2. Percentage contribution of climate variables under the Maxent modelling conducted in the present study. Definitions 
of variables are presented in Supplementary Materials, Table S3. Bold values refer to the most contributed variable in each 
species distribution modelling. 


Variable A. cruralis A. persica A. kermanensis C.agamuroides_ C.gastrophole  C. persepolense R. misonnei 
BIO2 5.9 8.7 15.7 

BIO4 8.5 39.4 44.5 

BIO6 34.3 21.3 51.2 16.4 46 
BIO9 19 55 

BIO11 23.6 28 40.6 9.8 9.6 

BIO12 20.3 5.2 19.9 42.9 
BIO13 12.2 15.7 

BIO14 6.6 5.4 

BIO15 13.8 6.1 

BIO16 16.6 

BIO17 15.9 6.4 

Slope 7.6 2.8 
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Figure 2. Bayesian Inference (Bl) gene tree of spider geckos inferred from 1670 bp of mitochondrial (12S and Cytb) and nuclear 
(MC1R) gene fragments. ML bootstrap support and posterior probability of Bayesian analyses are presented next to the nodes, 
respectively. Age estimated based on the substitution rates are denoted near the relevant nodes and include the mean and, 
between brackets, the HPD 95% confidence interval. 


Table 3. Mean + SD and range for significant characters of seven metric, meristic and ratio characters measured in Agamura 
population from Iranian Plateau. The right column refers to the significant values among populations. 


Character A. cruralis (n = 18) A. persica (n = 5) A. kermanensis (n = 7) P value 
Mean + SD Range Mean + SD Range Mean + SD Range 
HH 12.65+1.44 9.31-17.97 8.86+0.39 7.20-12.89 9.35+0.48 7.25-11.33 0.002 
SL 5.03+0.19 4.49-5.58 4.87+0.15 3.81-6.14 4.17+0.22 3.34-5.11 0.030 
10 10.11+0.41 8.81-11.40 8.79+0.29 6.60-11.22 10.85+0.68 8.59-13.63 0.005 
HLL 46.96+0.96 43.71-49.44 40.87+1.09 33.72-49.53 41.6643.36 30.03-56.67 0.042 
FLL/SVL 0.52+0.02 0.46-0.58 0.51+0.01 0.44-0.59 0.54+0.01 0.52-0.61 0.049 
NSA 28+0.94 25-30 37.11+1.01 31-47 27.1441.77 21-36 0.000 
NPV 53.20+2.10 46-59 54+0.77 49-59 48.85+1.07 46-54 0.014 
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Figure 3. Unrooted haplotype network of MC1R as nuclear marker. Circle size is proportional to the number of samples pre- 
haplotype, with colours corresponding to species in Fig. 2. Codes correspond to the specimens presented in Table S1. 


Table 4. Factor loadings of the first three principal 
components (PCs). The characters are defined in Table S4. 


Table 5. Factor loadings of the first three canonical variates 


(CVs). Characters were defined in Table S4. 


Characters PC1 PC2 PC3 Characters cv1 cv2 

HH 0.891 -0.079 0.172 FLLSVL -0.269 0.314 
SL 0.936 -0.005 0.034 NSA 0.701 -0.323 
10 0.857 0.298 0.110 NPV 0.110 0.209 
HLL 0.892 0.312 0.156 HH -0.320 0.980 
NSA -0.342 0.772 ~0.430 SL 0.922 -0.048 
FLLSVL -0.636 0.620 -0.056 10 -1.149 -1.103 
NPV 0.410 0.487 -0.748 HLL 0.492 0.706 
Eigenvalues 3.892 1.410 0.815 Eigenvalue 4.744 1.028 
Accumulated percent of trace 55.606 75.747 87.389 Accumulated percent of trace 82.2 100.0 


haplotypes were distinguished within the concatenated 
mitochondrial dataset and the nuclear MC1R marker 
included 32 unique haplotypes. 


Phylogenetic analyses and network construction 

The Bayesian and Maximum Likelihood trees showed 
identical topologies with high Bayesian posterior 
probabilities and bootstrap values (Figs. 2, $1). According 
to the phylogenetic analyses, Agamura is monophyletic 
(Fig. 2). Agamura is divided into three clades (Fig. 2) 
consisting of A. cruralis (the previously known species 
from eastern Iran), A. persica as a western clade (as 


mentioned by Szczerbak & Golubev, 1996) and A. 
kermanensis. Genetic distances (p-distance) of the two 
mitochondrial markers between clades reveals high 
diversity among the three lineages (i.e., 12S: 3.6-4.2%; 
Cytb: 8.3-13.2%). Variation within each clade is very low 
(12S: 0-0.8%; Cytb: 0.8-1.4%) especially within A. persica. 
The Cyrtopodion clade consists of four species that are 
distinctly separated from Agamura and each species is 
delimited from others with high bootstrap values and 
posterior probabilities support, but based on the Cytb 
tree (Supplementary Materials, Fig. S3) the genus was 
paraphyletic and Cyrtopodion gastrophole situated far 
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Figure 4. Potential species distribution models of spider geckos on the Iranian Plateau. A) C. agamuroides; B) A. persica; C) C. 
gastrophole; D) R. misonnei; E) A. cruralis; F) C. persepolense; G) A. kermanensis. The colours refer to the level of suitability as 
presented in the figure legends. 
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Figure 5. Ordination of principal component 1 (PC1) against PC2 (A) and canonical variate 1 (CV1) against CV2 (B) for all 
significant characters among Agamura clades in the Iranian Plateau 


from the genus Agamura. The Genus Rhinogekko is 
considered as sister genus to Bunopus and confirm the 
previous assumptions (Szczerbak & Golubev, 1996) (Fig. 
S3) and separated from the Agamura clade with high 
posterior probability and bootstrap support values. 

The haplotype network constructed for nuclear marker 
MC1R showed the similar pattern of a concatenated 
phylogenetic tree. According to the network (Fig. 3), 
separation of Rhinogekko and Cyrtopodion from the 
genus Agamura is confirmed, but the complexity in the 
genus Agamura remains. Agamura kermanensis is clearly 
differentiated from other clades with more mutations, 
but other two species have more similarity to each other. 


Divergence time estimates 
The analyses were run based on two mitochondrial genes 
(12S and Cytb) and finally, both trees linked together 
(Figs. 2, S2) and the results indicated that Rhinogekko 
split from angular-toed geckos around 17 million years 
ago (Mya; 95% HPD: 9.5-31.0 Mya). Divergence between 
Cyrtopodion and Agamura started through early-mid- 
Miocene ca. 15 Mya (95% HPD: 7.8-25.8 Mya). Speciation 
within the Cyrtopodion clade appears to have occurred 
11 Mya (95% HPD: 5.2-18.5 Mya). The split between C. 
gastrophole and C. agamuroides group took place around 
7 Mya (95% HPD: 3.5-12.7 Mya) and the separation 
between C. agamuroides and C. persepolense occurred 
2.8 Mya (95% HPD: 1.1-5.0 Mya). The cladogenesis of 
spider geckos started approximately in early Pliocene ca. 
4.3 Mya (95% HPD: 2.2-7.8 Mya) and the divergence of A. 
persica and A. cruralis have occurred 3.2 Mya (95% HPD: 
1.6-5.6 Mya). 

Spider geckos in the Iranian Plateau are distributed 
in all areas within the plateau across the mountains, 


plains and deserts. The ancestor of spider geckos was 
retrieved from the central part of the Iranian Plateau as 
the distribution of most of them covered the central part 
(Supplementary Materials, Fig. $4). Among spider geckos, 
Rhinogekko misonnei, Agamura persica, A. cruralis and 
A. kermanensis are distributed in the Central Plateau. 
Cyrtopodion species that were used in this study are from 
the Zagros area. It seems the spider gecko clade first 
appeared from the central part of the Plateau, because 
few of them are distributed in other parts. 


Species distribution modelling 

The results of Maxent modelling shows good AUCs for 
all models: A. cruralis AUC = 0.896+0.117; A. persica 
AUC = 0.879+0.028; A. kermanensis. AUC = 0.865+0.071; 
Cyrtopodion persepolense AUC = 0.997+0.001; C. 
agamuroides AUC = 0.864+0.188; C. gastrophole AUC = 
0.962+0.013; Rhinogekko misonnei AUC = 0.969+0.029. 
The predicted map indicated that suitable areas 
confirmed the current distribution of the species (Figs. 1, 
4). The most important climate variables are presented 
in Table 2. 


Morphological examination 

Twenty-seven morphological characters were examined 
among three distinct clades of the genus Agamura. 
Based on the analysis of variance (ANOVA) for metric 
characters, five characters were distinguished as 
significant characters (P < 0.05) (Head Height (HH), Snout 
Length (SL), Interorbital distance (10), Hind Limb Length 
(HLL) and Fore Limb Length/Snout-Vent length (FLL/ 
SVL)). The meristic characters were analysed by Kruskal- 
Wallis H test and two characters were distinguished 
as significant (Number of scales across widest part of 


abdomen (NSA), Number of scales between postmental 
scales and vent (NPV) (Table 3). PCA and CVA were 
calculated using significant characters obtained in the 
previous stage. In the PCA, the first three components 
explained 87.38% of total variance and in the CVA, the 
first two components explained 100% of total variance 
of characters (Fig. 5; Table 4, 5). 


DISCUSSION 


Phylogeny, diversity and endemism 

The molecular results confirmed the previously known 
and described species belong to the genus Cyrtopodion 
(Nazarov et al., 2009; Nazarov et al., 2012) and revealed 
three distinct lineages of the genus Agamura and 
confirmed their species status (Fig. 2). Agamura was 
defined as a monophyletic genus with four distinct 
species (Szczerbak & Golubev, 1996), but was recently 
revised morphologically and three of them were excluded 
from the genus (Anderson, 1999). Two subspecies were 
considered for Agamura persica: A. p. cruralis for the 
eastern population and A. p. persica for the western 
population, with some differences in morphological 
characters including the number of dorsal scales and 
shape of the tubercles (Anderson, 1999). Rhinogekko 
misonnei is a representative of the genus Rhinogekko 
and was distinctly separated from Agamura (high 
bootstrap support and posterior probability values). It 
has previously been considered in the genus Agamura by 
Szczerbak & Golubev (1996) and its distinction presented 
recently (Anderson, 1999; Krysko et al., 2007; Sindaco & 
Jereméenko, 2008). A high level of genetic differentiation 
between Rhinogekko and Agamura populations in our 
study strongly supports the exclusion of R. misonnei 
within the genus Agamura. Agamura gastropholis was 
one of the members of Agamura according to Szczerbak 
& Golubev (1996), but in the present study, this species 
was clustered within the genus Cyrtopodion. These 
geckos are common in the Iranian Plateau and have 
been recorded several times from arid regions of central 
and southern parts of Iran (Anderson, 1999; Nazarov 
et al., 2009; Moradi et al., 2011), but there were, until 
this study, no documented morphological or molecular 
evidence, and indeed no taxonomic study, on these 
geckos in the area. 

The Zagros Mountains are presented as an 
important region of endemism in west and southern 
Iran (Gholamifard, 2011; Hosseinzadeh et al., 2014). 
Recently, several species of lizards and snakes have been 
described, uncovering a deep history of radiations in the 
region (Nazarov et al., 2009; Ahmadzadeh et al., 2012; 
Rajabizadeh et al., 2012). All of the described species are 
endemic to the Iranian Plateau and their cladogenesis 
began in the Miocene. Interspecific variation of 12S and 
Cytb genes among the studied species is about 14-15% 
and 19-20% respectively, which is very high (Table 1). 
Species distribution modelling revealed that potential 
areas of distribution of these species did not have any 
overlap (Fig. 4), but different climate variables affect 
a species presence in a precise region (Hosseinian 
Yousefkhani et al., 2016). BIO6 (minimum temperature 
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in coldest month) is the important variable for the 
presence of three species (Cyrtopodion agamuroides, 
C. gastrophole and R. misonnei) (Table 2). Divergence 
between Cyrtopodion and Agamura took place 15 Mya 
by uplifting of the Zagros Mountains in the mid-Miocene 
and then cladogenesis within these clades occurred 
during the Pliocene (Fig. 2) with different climate 
conditions. 

Morphological variations between Agamura clades 
include the size of the body and especially in the length of 
limbs, which is longer in the eastern population than the 
western one. Ground structures in different parts of 
the range were observed directly during fieldwork and 
there are large rocks in the eastern and southern part 
of Iran, whereas the size of rocks decreases in central 
lran. Habitat features like larger rocks are related to the 
structure of adjacent mountains produced by geological 
events and pressures over millions of years (Anders et al., 
2010). 


Biogeography of spider geckos 

South-eastern and southern Iran are the probable 
regions for the spider gecko’s diversification. The collision 
of the Arabian and Eurasian Plates about 35 - 20 Mya 
(Dercourt et al., 1986; Mouthereau, 2011; McQuarrie & 
Van Hinsbergen, 2013) and the Zagros Mountain orogeny 
and consequent climate change are likely to have 
played an important role in diversification within spider 
geckos. According to our results, long branches with 
high bootstrap supports indicated a deep divergence of 
Rhinogekko from other spider geckos (17 Mya) and the 
Cyrtopodion from Agamura (15 Mya). But divergence 
time less than 4 Mya among A. persica and A. cruralis 
clades indicated to the historical taxon that originated 
in south-east Iran. Rhinogekko is restricted to the Lut 
Desert boundary and limited dispersal apparently caused 
by climate factors such as hot and dry conditions in the 
Lut area (Pourkhorsandi & Mirnejad, 2000). Closure 
of the new Tethys by collision with Arabia created the 
Zagros Basin (Khadivi, 2010), but contact with the Lut 
Block (Lut Block is a rigid plate on the Iranian Plateau) 
created a hot and dry region. The Cyrtopodion clade 
diversified in the late Miocene as a result of uplifting of 
the Zagros Mountains and aridification, in addition to 
trapping and isolating populations of this clade in the 
south Zagros valleys, which directly affected the variation 
among them. All species examined in this study are 
sensitive to temperature (Table 2) in different months 
of year except C. persepolense, which is dependent on 
annual precipitation. This isolation and diversification is 
directly reflected by the Zagros uplifting which affects 
precipitation in south and central Iran (Ramstein et al. 
1997). 

The main splits among spider geckos took place by 
vicariance in the mid-Miocene and other divergences 
within clades mainly represent different dispersal waves, 
aridification, climate change and restriction among 
valleys. Climate change and aridification directly affect 
R. misonnei; restriction among valleys is a common 
method for cladogenesis in the Cyrtopodion clade, as 
the newly described species belonging to this genus 
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are isolated among valleys and restricted to their small 
range. Cladogenesis in the genus Agamura may occurred 
by both vicariance events and dispersal waves on the 
Iranian Plateau from 15 Mya (Fig. 2). The southern clade 
(taxa. kermanensis) was isolated by the uplifting of the 
Lalezar Mountain in southern Kerman and the two others 
(A. persica and A. cruralis) diverged from each other 
about 3 Mya by occupying different niches. 


Taxonomic implications 

The monophyly of Cyrtopodion and Agamura is 
confirmed and Rhinogekko is placed at the root of the 
tree as the first diverged spider geckos in south-east 
Iran. However, the Agamura divided into three clades 
that are phylogenetically distinct from others, and each 
clade might represent a distinct species (Fig. 2). This is 
confirmed by our analyses of the three major clades of 
Agamura represented on the Iranian Plateau; two of 
them refer to the species previously known as A. cruralis 
(eastern clade, including the type locality of A. cruralis 
from southern Baluchestan, Bahukalat) and A. persica 
(western clade) (Hora, 1926; Szczerbak & Golubev, 1996). 
The third clade in Agamura (A. kermanensis) diverged 
widely from the two others (13.2% and 11.7%; genetic 
differentiation from A. cruralis and A. persica in Cytb, 
respectively). Agamura kermanensis can be considered 
as a new taxon based on morphological and molecular 
evidence and will be described soon by authors. The 
morphological variation between the newly explored 
clade and two other Agamura clades is clearly explained 
in Figure 5. Species niche modelling confirmed this 
separation by estimating separate suitable areas (Fig. 
4). The haplotype network based on the nuclear marker 
clearly showed separation of Agamura kermanensis (Fig. 
3). 


CONCLUSIONS 


The present study provides a biogeographical view of the 
spider geckos in the Iranian Plateau. The phylogenetic 
history and divergence times of spider geckos support 
the geological events during the Miocene period. The 
monophyly of Agamura and Cyrtopodion indicate an 
important role of Zagros orogeny in the isolation of these 
lineages during the mid-Miocene. Following dispersal, 
Agamura rapidly diversified in the Iranian Plateau. The 
southern populationin Kerman province showed a distinct 
clade from other Agamura populations. The southern 
Zagros slopes were involved in different collision events 
with Arabia and then new valleys appeared. Radiation 
within Cyrtopodion clade refers to these valley creations 
taking place during the late Miocene. Rhinogekko is a 
distinct clade that first diverged from all other spider 
geckos by closing the new Tethys about 17 Mya in the 
Lut region. Finally, our phylogenetic analyses suggest 
that morphological diversity among these geckos arose 
by historical process such as Zagros orogeny, which 
can explain the basal divergence within the group. The 
southern clade will be described as a new taxon ina 
separate article by authors in future. 
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